Available data on Makay protected area

Author

Florent Bédécarrats

Introduction

This interactive workbook aims at providing an overview of data sources that are potentially relevant for the analysis of economic, social and ecological dynamics in the area concerned by the creation of the Makay protected area.

Working environment

We use R and a series of packages dedicated to spatial data and analysis.

Code
library(aws.s3)
library(tidyverse) # toolkit for data manipulation
library(geodata) # to get municipalities
library(sf) # to handle spatial data (vectors)
library(terra) # to handle patial data (rasters)
library(mapme.biodiversity) # to compute spatial indicators
library(tmap) # nice maps
library(zoo) # time series
library(units) # ensures units are right
library(future) # to parallelize computations
library(exactextractr) # engine for mapme.biodiversity
library(SPEI) # to compute rainfall  

Makay protected area

We load the Makay protected area boundaries.

Code
makay_ap <- st_read("data/Makay_wgs84.geojson", quiet = TRUE) %>%
  rename(Statut = SOURCETHM) %>%
  mutate(Statut = recode(Statut, "zone tampon" = "Zone tampon"))

# Fist we dissolve the different areas
makay_union <- makay_ap %>%
  st_union() %>%
  st_sf() %>%
  st_make_valid()

# Draw map
tmap_mode("view")
tm_shape(makay_ap) +
  tm_polygons(col = "Statut", alpha = 0.5) +
  tm_shape(makay_union) +
  tm_borders(col = "darkblue") + 
  tm_add_legend(type = "fill", col = "darkblue", labels = "Limite externe") +
  tm_basemap("OpenStreetMap")

Makay protected area (source: Auréa?)

The protected area is composted of five polygons with different status: four polygons are cores, and one polygon is a buffer area.

Administrative areas

We will now select all municipalities located within a 20km distance from the outer border of the Makay protected area (core or buffer zone).

Code
buff_size <- 20
# Create a 20km buffer around Makay PA
makay_ext_buffer <- makay_union %>%
  st_buffer(dist = set_units(buff_size, km)) %>%
  st_as_sf()
# Download Madagascar municipalities
muni_mada <- gadm("MDG", level  =4, path = "data") %>%
  st_as_sf() %>%
  rename(Nom_province = NAME_1, Nom_region = NAME_2, Nom_district = NAME_3,
         Nom_commune = NAME_4)
# Select municipalities that intersect the buffer
muni_makay <- muni_mada %>%
  st_filter(makay_union, .predicate = st_intersects)
muni_makay_names <- pluck(muni_makay, "Nom_commune")
muni_around_makay <- muni_mada %>%
  st_filter(makay_ext_buffer, .predicate = st_intersects)
muni_around_makay_names <- pluck(muni_around_makay, "Nom_commune")
muni_makay <- muni_mada %>%
  st_filter(makay_union, .predicate = st_intersects)

muni_around_ext_makay <- muni_around_makay_names[!(muni_around_makay_names %in% 
                                                     muni_makay_names)]

# Draw map
tm_shape(makay_union) +
  tm_borders(col = "darkblue") + 
  tm_shape(muni_around_makay) + 
  tm_borders(col = "purple") +
  tm_fill(alpha = 0, id = "Nom_commune", 
          popup.vars = c("Nom_commune","Nom_district", "Nom_region", 
                         "Nom_province")) +
  tm_basemap("OpenStreetMap")

Surrounding municipalities (source: GADM)

There are 13 municipalities that directly overlap with Makay: Tandrano, Beroroha, Fanjakana, Mandronarivo, Marerano, Sakena, Tanamary, Ambia, Ankilizato, Beronono, Malaimbandy, Mandabe, Tsimazava. There are 4 municipalities that do not directly overlap with Makay PA but are located within a 20 km radius from Makay PA: Mandrosonoro, Fitampito, Behisatra, Beharona.

Rainfall

We use the CHIRPS data from NASA to estimate 3 day average rainfall on the Makay protected area.The source data has a spatial resolution of 0.05° (about 5km). We compute the data using the package mapme.biodiversity (Görgen and Bhandari 2023).

Code
# We group Makay PA zones and communes to have an object with all areas of interest
# We keep only one common field "name"
aoi <- bind_rows(
  filter(makay_ap, Statut == "Zone tampon"),
  filter(makay_ap, Statut == "Noyau dur") %>% 
    mutate(Statut = paste0(Statut, "_", 1:4))) %>%
  select(name = Statut) %>%
  bind_rows(
    mutate(makay_union, name = "Ensemble")) %>%
  mutate(name = paste0("AP Makay : ", name)) %>%
  bind_rows(select(muni_around_makay, name = Nom_commune)) %>%
  st_cast("POLYGON")

# we use parallel computing to reduce processing time
plan(multisession, workers = 8)

# Reference portfolio (chirps data was loaded at startup)
aoi <- init_portfolio(aoi, years = 1990:2020,
                      outdir = "data") 


aoi_1 <- s3read_using(FUN = read_rds, object = "diffusion/makay/aoi_1.rds",
                      bucket = "fbedecarrats", opts = list("region" = ""))

if (!exists("aoi_1")) {
  aoi <- aoi %>%
  get_resources("chirps")
  # Compute precipitation
  aoi_1 <- aoi %>%
    calc_indicators("precipitation_chirps",
                    engine = "exactextract",
                    scales_spi = 3,
                    spi_prev_years = 8)
  s3write_using(aoi_1, FUN = write_rds, object = "diffusion/makay/aoi_1.rds",
                bucket = "fbedecarrats", opts = list("region" = ""))
}

# We compute 3d average precipitations
precip_3d <- aoi_1 %>%
  unnest(precipitation_chirps) %>%
  group_by(name) %>%
  mutate(rainfall_3d_mean = rollmean(absolute, k = 3, fill = NA))

# We display precipitations for Makay area
precip_3d %>%
  filter(name == "AP Makay : Ensemble") %>%
  ggplot(aes(x = dates, y = rainfall_3d_mean)) +
  geom_line()

Rainfalls on Makay PA (source : CHIRPS, 3 day rolling average)

Cyclones

We use the Tropical cyclone best track data (IBTrACS) to get history of cyclone trajectories and cyclone attributes (Knapp et al. 2010), as well as raw code from (Schielein 2023).

Code
if (!file.exists("data/ibtracs/IBTrACS.since1980.list.v04r00.lines.shp")) {
  dir.create("data/ibtracs")
  download.file(url = "https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/shapefile/IBTrACS.since1980.list.v04r00.lines.zip", 
                destfile = "data/ibtracs/ibtracs_lines.zip")
  unzip("data/ibtracs/ibtracs_lines.zip", exdir = "data/ibtracs")
}

mada <- gadm("MDG", level=0, path = "data") %>%
  st_as_sf()

# column description is given here: https://www.ncei.noaa.gov/sites/default/files/2021-07/IBTrACS_v04_column_documentation.pdf
cyclones <- read_sf("data/ibtracs/IBTrACS.since1980.list.v04r00.lines.shp")

# Create a buffer 300km around Madagascar
mada_buff <- st_buffer(mada, 300)

cyclones$wind_combined <- cyclones %>%
  select(contains("WIND")) %>%
  select(-WMO_WIND) %>%
  st_drop_geometry() %>%
  rowMeans(., na.rm = T)

cyclones_mada <- cyclones %>%
  st_intersection(mada_buff) %>%
  select(Season = SEASON, Name = NAME, Winds_in_knots = wind_combined)

tm_shape(makay_union) + 
  tm_borders(col = "darkblue") + 
  tm_shape(cyclones_mada) +
  tm_lines(col = "Season", lwd = "Winds_in_knots", scale = 3, palette = "Blues",
           popup.vars = c("Name", "Season", "Winds in knots" = "Winds_in_knots"),
           legend.format = list(big.mark=""), popup.format = list(big.mark=""))

Cyclone trajectories, dates and intensity (source: IBTrACS)

The map shows that the Makay protected area has been hit by a number of significant cyclones. Most notably, Freddy in 2023 and Bastiraï in 2022, both with winds of around 57 knots (1 knot = 1.825 km/h). Lesser cyclones had hit the area in previous years: Chedza in 2014 (winds around 42 knots), Elita in 2005 (winds around 52 knots).

Fires

This analysis is based on the example produced by Darius Goergen on the Mapme website.

Code
aoi_2 <- get_resources(aoi_1, "nasa_firms", instrument = "MODIS")

gpkgs <- list.files(file.path("data", "nasa_firms"),
                    pattern = ".gpkg", full.names = TRUE)
nasa_firms <- map_dfr(gpkgs, function(x) {
  read_sf(x, wkt_filter = st_as_text(st_as_sfc(st_bbox(aoi_2))))
})
nasa_firms <- nasa_firms[unlist(st_contains(aoi_2, nasa_firms)), ]
nasa_firms <- filter(nasa_firms, confidence > 50)
nasa_firms$year <- as.factor(year(nasa_firms$acq_date))
# Remove year 2000 for visibility
nasa_firms <- filter(nasa_firms, year != "2000")

tmap_mode("plot")

tm_shape(nasa_firms) +
  tm_dots(col = "darkorange", alpha = 0.1) + 
  tm_facets(by = "year", nrow = 5, ncol = 4) +
  tm_shape(makay_union) + 
  tm_borders(col = "darkblue")

Fires in and around Makay PA (source: MODIS)

Compute indicators on fires.

Code
aoi_2 <- aoi_2 %>%
  calc_indicators("active_fire_counts")

fires <- aoi_2 %>%
  unnest(active_fire_counts)

fires %>%
  filter(name == "AP Makay : Ensemble" & year != 2000) %>%
  ggplot(aes(x = year, y = active_fire_counts, group = name)) + 
  geom_line(col = "darkorange")

Fire occurence count in Makay PA (source: MODIS)

Forest cover and loss

This data relies on global forest cover.

Code
aoi_2 <- aoi_2 %>%
  get_resources(c("gfw_treecover", "gfw_lossyear"))
tc_2000 <- rast("data/gfw_treecover/Hansen_GFC-2021-v1.9_treecover2000_20S_040E.tif") %>%
  crop(aoi_2)
  

tm_shape(tc_2000) +
  tm_raster(title = "Forest cover in 2000 (%)", palette = "YlGn", style = "cont") +
  tm_shape(makay_union) + 
  tm_borders(col = "darkblue") + 
  tm_layout(legend.outside = TRUE) + 
  tm_add_legend(type = "line", labels = "Makay PA delimitation", col = "darkblue")

Forest cover in 2000 (source: GFW)

Forest cover loss

Code
lossyear <- rast("data/gfw_lossyear/Hansen_GFC-2021-v1.9_lossyear_20S_040E.tif") %>%
  crop(aoi_2)
loss <- lossyear
values(loss) <- ifelse(values(lossyear) > 0, 1, NA)

tm_shape(tc_2000) +
  tm_raster(title = "Forest cover in 2000 (%)", palette = "YlGn", style = "cont") +
  tm_shape(loss) + 
  tm_raster(title = "Forest cover loss 2001-2021", 
            style = "cat", breaks = c(0, 0.5, 1.5), palette = "red") +
  tm_shape(makay_union) + 
  tm_borders(col = "darkblue") + 
  tm_layout(legend.outside = TRUE) + 
  tm_add_legend(type = "line", labels = "Makay PA delimitation", 
                col = "darkblue")

Forest cover in 2000 and loss since (source: GFW)

Statistincs on forest cover loss

Code
if (Sys.getenv("USER") == "onyxia") {
  aoi_3 <- s3read_using(FUN = read_rds, object = "diffusion/makay/aoi_3.rds",
                      bucket = "fbedecarrats", opts = list("region" = ""))
} else {
  aoi_3 <- aoi_2 %>%
    calc_indicators("treecover_area", min_cover = 10, min_size = 0.5)
  s3write_using(aoi_3, FUN = write_rds, object = "diffusion/makay/aoi_3.rds",
                bucket = "fbedecarrats", opts = list("region" = ""))
}

treecover <- aoi_3 %>%
  unnest(treecover_area)

treecover_makay <- treecover %>%
  filter(name == "AP Makay : Ensemble")

treecover_makay_2000 <- treecover_makay %>%
  filter(years == 2000) %>%
  pluck("treecover") %>%
  round()
treecover_makay_2020 <- treecover_makay %>%
  filter(years == 2020) %>%
  pluck("treecover") %>%
  round()
varcover <- (treecover_makay_2020 - treecover_makay_2000) / treecover_makay_2000

treecover_makay %>%
  ggplot(aes(x = years, y = treecover, group = name)) + 
  geom_line(col = "darkgreen")

Note that the amont loss looks great because the graph origin is not 0. In reality the observed loss goes from hectares in 2000 to hectares in 2020, that is a variation of -4.88 % in 20 years.

With an origin at 0, here is what the graph looks like:

Code
treecover_makay %>%
  ggplot(aes(x = years, y = treecover, group = name, ymin = 0)) + 
  geom_line(col = "darkgreen")

Droughts

Conflicts

ACLED data

We use here the data for all type of conflicts from ACLED (Armed Conflict Location & Event Data Project (ACLED); www.acleddata.com), for which access was requested on August 17th, 2023. ACLED data records violent events registered between 1997 and today. For Madagascar, the dataset includes 2,079 events (761 battles, 639 riots, 50 explosions/remote violence and 629 events of violence against civilians), corresponding to 3,418 fatalities. There seems to be a strong biais towards recent events:

Events registered in Acled for Madagascar per year

Click here to explore on ACLED website the content of their database for Madagascar since 1997.

UCDP data (1989-2022)

More precisely, this is the UCDP Georeferenced Event Dataset (GED) Global version 23.1

“This dataset is UCDP’s most disaggregated dataset, covering individual events of organized violence (phenomena of lethal violence occurring at a given time and place). These events are sufficiently fine-grained to be geo-coded down to the level of individual villages, with temporal durations disaggregated to single, individual days.”

Code
ucd <- s3read_using(FUN = read_rds, 
                    object = "diffusion/makay/GEDEvent_v23_1.rds",
                    bucket = "fbedecarrats", opts = list("region" = "")) %>%
  filter(country == "Madagascar (Malagasy)") %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

tmap_mode("view")
tm_shape(ucd) +
  tm_dots()

It appears that there are very few occurrences referenced in this dataset.

There are 47 occurrences registered in UCDP dataset. It appears to be very much focused on violent event involving public actors (“organized” violence).

Species occurrences

Use of the GBIF data. GBIF.org (17 August 2023) GBIF Occurrence Download https://doi.org/10.15468/dl.w6mtcs

Code
# Waiting for data access

References

Görgen, Darius A., and Om Prakash Bhandari. 2023. “Mapme.biodiversity: Efficient Monitoring of Global Biodiversity Portfolios.”
Knapp, Kenneth R., Michael C. Kruk, David H. Levinson, Howard J. Diamond, and Charles J. Neumann. 2010. “The International Best Track Archive for Climate Stewardship (IBTrACS) Unifying Tropical Cyclone Data.” Bulletin of the American Meteorological Society 91 (3): 363–76.
Schielein, Johannes. 2023. “Mapme.protectedareas: Reproducible Processing and Analysis Routines with r to Assist in Planning Monitoring and Evaluation of Projects to Support PAs Worldwide.”